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ABSTRACT 


The  objective  of  this  study  was  to  develop  an 
analytical  model  of  the  basic  hierarchical  structure  of  the 
human  bone.  The  model  computed  the  stiffness  of  composite 
collagen  fibers  comprised  of  collagen  fibrils  and 
hydroxyapatite  mineral  crystals.  Next,  the  stiffness  of  the 
concentric  lamella  was  computed  utilizing  the  stiffness  of 
the  collagen  fibers  and  layer  information.  Finally,  the 
effective  stiffness  of  the  bone  was  estimated.  In  order  to 
determine  the  stiffness  of  the  collagen  fiber,  a  three- 
dimensional  finite  element  model  was  developed  and  a  simple 
analytical  model  was  derived.  The  simple  analytical  model 
was  validated  using  the  finite  element  results.  The 
lamination  theory  of  unidirectional  fibrous  composites  was 
used  to  calculate  the  stiffness  of  the  lamella  and 
eventually  the  bone  stiffness.  A  series  of  parametric 
studies  were  conducted  to  understand  what  parameter (s) 
affected  the  stiffness  of  the  bone  most  significantly. 
This  information  will  be  useful  when  an  artificial  bone 
structure  is  designed. 
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I .  INTRODUCTION 


This  thesis  research  developed  both  an  analytical 


model  and  a 

finite 

element 

model 

of  the 

hierarchical 

structure  of 

human 

bone . 

The 

analytical 

model  was 

validated  using  the  results  from  the  finite  element  model. 
Before  we  delve  into  the  aspects  of  each  model,  it  is 
important  to  understand  the  structure  of  the  human  bone. 
Figure  1  shows  molecular  to  ultra  to  micro  to  macro 
components  of  the  bone,  which  will  be  discussed  in  the 
subsequent  section. 


Figure  1.  Hierarchical  Structure  of  Human  Bone  (From: 

Rockwood,  2001) 
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A.  HUMAN  BONE  STRUCTURE 

Starting  with  the  molecular  structure  of  bone  we  have 
hydroxyapatite  mineral  crystals,  commonly  known  as  apatite. 
These  mineral  crystals  provide  for  the  hardness  and  brittle 
nature  of  bone  (Apatite,  2001)  . 

The  hydroxyapatite  lies  in  sections  around  a  collagen 
fibril.  The  collagen  has  great  tensile  strength  and  allows 
for  bone's  elastic  behavior  which  complements  the 
properties  of  the  apatite  (Collagen,  2001) .  Together  these 
materials  make  up  the  composite  fiber  and  account  for  the 
molecular  structure  of  bone. 

The  ultra  component  of  the  bone  structure  consists  of 
many  concentric  lamellas.  In  nature  these  lamella  are 
oriented  45°  from  the  longitudinal  axis.  Each  layer  of  the 
lamellae  are  at  opposing  angles  and  there  are  3-20  lamella 
comprising  the  lamellae  (Caceci,  2001). 

The  osteon  or  the  Haversion  system  is  the  micro 
component  of  bone.  This  is  made  up  of  many  lamellae 
encircling  an  inner  region  to  provide  protection.  This 
inner  region  consists  of: 

-  The  Haversion  canal  which  contains  one  or  two  small 
blood  vessels  and  a  nerve  (Osteon,  2002) . 

Lacunae;  which  are  holes  for  one  bone  cell 
(osteocyte)  to  live  (Osteon,  2002) . 

Canaliculi  which  enable  osteocytes  to  communicate 
with  one  another  (Osteon,  2002)  . 

-  Volkman' s  canal  which  enables  adjacent  osteons  to 
communicate  with  each  other  (Osteon,  2002) . 
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The  macro  components  of  bone  make  up  bone  as  most  of 
us  know  it  today.  As  seen  in  Figure  1,  there  are  many 

different  parts  which  make  up  the  outer  structure  of  the 
bone.  This  research  did  not  account  for  every  one,  but 

instead  grouped  them  together  to  account  for  the  30%  of 
bone  which  was  not  analyzed  with  the  finite  element  model 
or  with  the  analytic  predictions.  In  other  words  our  final 
values  for  the  properties  of  bone  were  based  upon  the 

lamella  structure  and  multiplied  by  70%. 

B.  OBJECTIVES 

The  objective  of  this  study  was  to  develop  an 
analytical  model  of  the  basic  hierarchical  structure  of  the 
human  bone,  consisting  of  collagen  and  hydroxyapatite.  The 
model  computed  the  stiffness  of  composite  collagen  fibers 

comprised  of  collagen  fibrils  and  hydroxyapatite  mineral 
crystals.  Once  the  model  was  developed,  the  stiffness  of 
the  concentric  lamella  was  computed  utilizing  the  stiffness 
of  the  collagen  fibers  and  layer  information.  Finally,  the 
effective  stiffness  of  the  bone  was  estimated. 

In  order  to  determine  the  stiffness  of  the  collagen 
fiber,  a  three-dimensional  finite  element  model  was 
developed  and  a  simple  analytical  model  was  derived.  The 
simple  analytical  model  was  validated  using  the  finite 
element  results.  The  lamination  theory  of  unidirectional 
fibrous  composites  was  used  to  calculate  the  stiffness  of 
the  lamella  and  eventually  the  bone  stiffness.  A  series  of 
parametric  studies  were  conducted  to  understand  what 
parameter (s)  affected  the  stiffness  of  the  bone  most 
significantly . 
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II.  FINITE  ELEMENT  MODEL 


A.  MODEL  OVERVIEW 

A  three-dimensional,  fixed  displacement,  finite 
element  model  of  the  most  basic  human  bone  structure  was 
developed  that  is  capable  of  predicting  the  total  stress  of 
the  composite  fiber  and  therefore,  through  analytical 
calculat ions ,  the  elastic  modulus  of  the  composite  fiber. 
The  model  is  on  a  nanometer  scale  and  represents  one 
section  of  the  composite  fiber,  including  collagen  and 
hydroxyapatite  mineral  crystals.  The  commercial  finite 
element  package  MSC/PATRAN  was  used  for  pre  and  post 
processing  and  MSC/NASTRAN  was  used  for  analysis. 

B .  GEOMETRY 


One 

section  of 

the 

composite  fiber  was  modeled. 

as 

shown  in 

Figure  2 . 

Beginning  with  a 

solid  cylinder 

to 

represent 

a  portion 

of 

the  collagen 

fibril,  64  nm 

in 

length  with  a  radius  of  60  nm,  hydroxyapatite  mineral 
crystals,  48  nm  in  length  and  2nm  in  width,  were  placed  on 
the  outside  of  the  cylinder.  These  two  materials  made  up 
the  composite  collagen  fiber.  The  elastic  modulus  of  each 
component  was  varied  in  order  to  analyze  the  individual 
effect  on  the  elastic  modulus  of  the  composite.  The 
stiffness  of  collagen  was  adjusted  from  10  to  30  GPa  at  5 
GPa  intervals.  Separately,  the  stiffness  of  hydroxyapatite 
was  adjusted  from  700  GPa  to  450  GPa  at  50  GPa  intervals. 
These  values  were  chosen  at  random  to  determine  the  trend 
of  the  composite  stiffness  as  each  value  was  changed. 
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Figure  2 . 


Finite  Element  Model  of  Composite  Fiber 


C .  STATIC  RESPONSE 

The  finite  element  model  was  run  multiple  times  with 
different  values  for  the  elastic  modulus  of  collagen  and 
hydroxyapatite  to  determine  which  material  affected  the 
elastic  modulus  of  the  composite  most  dramatically. 
Applying  boundary  conditions  enabled  us  to  fix  the 
translation  in  the  x-direction,  fix  the  translation  in  the 
y-direction,  and  fix  rotation  in  all  directions.  Applying 
a  fixed  displacement  to  the  model  in  the  z-direction,  we 
were  able  to  determine  the  stress  at  each  point  of  the 
composite  fiber.  Through  analytic  calculations,  we  then 
determined  the  average  stress  over  the  cross  sectional 
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area 


to  ultimately  obtain  the  elastic  modulus  of  the 


composite  fiber 


where  E 


elastic  modulus,  a  = 


stress,  and  s  =  strain.  Once  these  results  were  tabulated 
we  used  them  to  validate  our  analytic  model. 


Figure  3  is  an  example  of  the  output  of  the  model  when 
Eapatite  =  450  GPa  and  Ecoliagen  =  20  GPa.  The  various  colors 
represent  the  stresses  on  the  given  portion  of  the 
composite  material.  The  analytic  calculations  to  determine 
the  average  stress  were  accomplished  on  the  cross  sectional 
area  (face  of  the  fiber)  .  The  model  is  on  a  nanometer 
scale,  so  to  obtain  a  meaningful  value,  we  need  to  divide 
each  stress  value  by  the  average  stress  over  the  cross 
sectional  area  to  obtain  the  stress  concentration.  For 
example,  the  white  value  for  stress  is  2.79e-ll  N/nm2  which 
when  converted  to  common  dimensions,  becomes  2.79e7  N/m2  or 
27.9  MPa.  We  divide  this  value  of  stress  by  the  average 
stress  obtained  through  earlier  calculations  and  obtain, 
6.50  stress  per  unit  stress,  the  stress  concentration. 
This  value  is  explained  to  give  the  reader  a  general  idea 
of  what  they  are  seeing  in  Figure  3,  and  will  make  more 
sense  when  we  discuss  our  results. 
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Figure  3. 


Stress  over  Composite 


From  the  finite  element  model,  we  gained  insight  into 
the  composite  fiber  properties  and  used  the  results  to 
validate  a  simple  analytic  model. 
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Ill . 


ANALYTIC  SOLUTION 


An  analytic  model  was  derived  in  a  sequence  of  three 
separate  computer  programs.  The  first  computed  the 
stiffness  of  composite  collagen  fibers  comprised  of 
collagen  fibrils  and  hydroxyapatite  mineral  crystals.  Once 
this  analytic  model  was  validated  against  the  results  of 
the  finite  element  model,  the  stiffness  of  the  concentric 
lamella  was  computed  utilizing  the  stiffness  of  the 
collagen  fibers  and  layer  information.  Finally,  the 
effective  stiffness  of  the  bone  was  estimated. 

A.  FORCE-DEFORMATION  BEHAVIOR  OF  A  TYPICAL  UNIFORM  AXIAL- 

DEFORMATION  ELEMENT 

An  analytic  model  was  derived  and  the  results  from  the 
finite  element  model  were  used  as  validation  for  this 
model.  As  a  starting  point,  a  MATLAB  computer  program 
using  uniform  axial  deformation  (spring  constant)  theory 
was  developed.  It  enabled  the  user  to  change  the  values 
for  the  elastic  modulus  of  hydroxyapatite  and  collagen  to 
determine  the  longitudinal  stiffness  of  the  composite 
fiber.  Figure  4  shows  a  schematic  of  the  composite  fiber. 
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Figure  4.  Spring  Constant  Schematic  of  Composite  Fiber 
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We  can  compare  Figure  4  to  Figure  5  to  see  how  the 


spring  model  correlates  to  the  composite  model. 


SECTION  1  SECTION  2  SECTION  3 


1-2  L  3 

Figure  5.  2-D  Finite  Element  Model  Schematic  Showing 

Sections 

The  basic  theory  behind  this  solution  is 

k  kt  k2  k2 

termed  the  force-deformation  behavior  of  a  typical  uniform 
axial-deformation  element  (Craig,  2000);  where  k  =  the 
total  stiffness  coefficient  (force  required  to  produce  a 
unit  elongation  of  the  member)  ,  ki  =  the  stiffness 
coefficient  for  the  first  section  of  composite,  k2  =  the 
stiffness  coefficient  for  the  second  section  of  the 
composite,  and  k3  =  the  stiffness  coefficient  for  the  third 
section  of  the  composite.  Note  that  k2  =  k2a  +  k2c  and  that 
k2a  =  hydroxyapatite  stiffness  coefficient  for  the  second 
section  and  k2c  =  the  collagen  stiffness  coefficient  for  the 
second  section.  Solving  this  equation  for  k  we  find 
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this,  we  find  that  =  (a  +  dcr)s7  -  as2 -rsdz  =  0.  Simplifying 

T  z, 

this  equation  and  integrating  we  find  that  cr(z)  =  cro+—  where 
cro=0  as  there  is  no  stress  at  the  end  of  the  cross  section. 


tz 

Assuming  that  r  is  constant,  cr  =  — 

s 


occurs  at  the  center  of  the  section. 


and  the  maximum  stress 
as  seen  in  figure  7, 
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Figure  7 . 


Stress  Over  One  Section  of  Hydroxyapatite 
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Figure  8.  Shear  Stress  Over  One  Section  of 

Hydroxyapatite 
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From  figures  7  and  8  we  find  that  Fnv=crnvA  = —  =  J ,eai 
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From 


this  equation  we  determine  that  our  value  for  k2a  must  have 
a  factor  of  0.25  in  the  equation,  therefore, 
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where  r  =  radius,  the  c  subscript  represents  collagen,  the 
a  subscript  represents  apatite,  and  n  =  the  number  of 
sections  of  apatite  around  the  fibril.  We  can  now 
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determine  k,  which  in  turn  enables  us  to  determine  the 
elastic  modulus  of  the  composite  fiber,  using 

_  kL 

^ composite  2 

nr 

The  value  obtained  here  for  ECOmposite/-  is  the 
longitudinal  elastic  modulus  of  the  composite  fiber.  The 
results  from  this  analytic  process  were  used  to  validate 
our  model  and  as  a  basis  to  move  on  to  the  next  step  which 
was  to  find  the  elastic  modulus  of  the  composite  fiber  when 
oriented  at  various  angles. 


B.  LAMINATION  THEORY  OF  UNIDIRECTIONAL  FIBROUS  COMPOSITES 

The  lamination  theory  of  unidirectional  fibrous 
composites  was  utilized  to  derive  the  next  portion  of  the 
analytic  model.  The  value  for  the  longitudinal  elastic 
modulus  of  the  composite  fiber  was  the  starting  point. 
This  value  was  used  in  the  transformed  stiffness  matrix  to 
produce  the  elastic  modulus  of  a  single  fiber  oriented  at 
various  angles.  The  angles  of  orientation  included  15°, 
30°,  45°,  60°,  -15°,  -30°,  -45°,  and  -60°. 

The  following  stress-strain  relations  for  isotropic 
materials  were  used  (Jones,  1975)  : 


in  which  Q{j  are  the  reduced  stiffnesses;  specifically, 

E 

Qn= - ^  where  v  is  Poisson's  ratio  and  assumed  to  be  0.3, 

1  -v 

and  Qn ,  Q22  ,  and  Q66  are  all  zero.  Using  these  results,  we 
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then  calculated  the  transformed  reduced  stiffnesses  using 
the  following  stress-strain  relations  in  the  xy  coordinate 
directions  (Jones,  1975): 


Qu 

Ql2 

Qi«' 

.  = 

Q\2 

Q22 

Q26 

< 

.V 

_Qi« 

Q26 

Qbb  _ 

A. 

in  which  Qtj  are  the  transformed  reduced  stiffnesses; 
specifically  (Jones,  1975), 

Qu  =  Qn  cos4  0  +  2 [Qn  +  2 Qb6 ) sin2  0 cos2  0  +  Q22  sin4  0 
Qn  =  {Qn  +  Q22  ~  4266 ) sin2  e c°s2  0  +  QX2  (sin4  0  +  cos4  0 ) 

Q22  =  Qi 1  sin4  0  +  2 (Ql2  +  2Q66 )  sin2  0 cos2  0  +  Q22  cos4  0 

Qn  =  (2i  1  -  fin  -  2266 )  sin  0 cos3  0  +  (Ql2  -  Q22  +  2 Q66 )  sin3  0  cos  0 

Q26={Qn- Qn  ~ 2Q66) sin3  6* cos 0  +  ( Ql2  - Q22  +  2 Qbb ) sin 0 cos3  0 

Gee  =  (Qn  +  Q22  ~  2Qn  -  2Qbb ) sin2  0 cos2  0  +  Q66  (sin4  0  +  cos4  0 ) 

The  transformed  elastic  modulus  of  the  composite  fiber 
was  equivalent  to  Qn  .  This  value  was  then  used  as  a 
starting  point  in  the  last  step  of  the  analytic  model 
derivation  to  determine  the  actual  elastic  modulus  of  human 
bone . 

C.  LAMELLAE  AND  BONE  ELASTIC  MODULUS 

Considering  actual  human  bone  structure,  the  composite 
fibers  were  wrapped  in  a  concentric  lamella,  each  oriented 
at  opposing  angles.  From  the  previous  section,  we 

determined  the  elastic  modulus  of  a  given  composite  fiber 
oriented  at  a  given  angle.  The  last  step  in  the  sequence 
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to  derive  the  analytic  model  started  with  these  results  and 
determined  first,  the  elastic  modulus  of  the  lamellae 
structure,  and  second,  the  elastic  modulus  of  bone. 

The  lamellae  was  assumed  to  have  12  layers,  as  an 
average  between  the  3-20  lamella  stated  in  the  references 
(Caceci,  2001)  .  The  transformed  elastic  modulus  for  each 
angle  (15°,  30°,  45°,  60°,  -15°,  -30°,  -45°,  and  -60°),  from 
the  previous  derivation  using  the  transformed  matrix,  was 

the  starting  point.  Radius  and  area  (A  =  nr~^ calculations 

were  made  for  each  layer  of  the  concentric  lamella  as  well 
as  for  the  total  area  of  the  lamellae.  To  calculate  the 
elastic  modulus  of  the  entire  lamellae,  the  following 
formula  was  used: 


(=i 


TotalArea 

Once  the  stiffness  of  the  lamellae  was  obtained,  the  only 
thing  which  remained  was  to  make  an  estimate  to  determine 
the  actual  stiffness  of  the  entire  bone.  To  do  this,  we 
assumed  that  the  lamellae  represented  70%  of  the  human  bone 
structure.  We  multiplied  the  lamellae  stiffness  by  0.7  to 
account  for  this  70%  and  obtained  the  elastic  modulus  of 
human  bone . 
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IV.  RESULTS  AND  DISCUSSION 


A.  EFFECT  OF  ELASTIC  MODULUS  OF  COLLAGEN  AND 

HYDROXYAPATITE  ON  STRESS  CONCENTRATION  OF  COMPOSITE 
FIBER 

From  the  figures  below  we  determine  that  the  elastic 

modulus  of  collagen  has  a  similar  impact  on  the  stress 

concentration  of  the  composite  fiber  as  the  elastic  modulus 
of  hydroxyapatite,  although  collagen  may  be  somewhat  more 

influential.  The  stress  results  of  the  finite  element 
model  enabled  us  to  determine  the  stress  concentration  over 
the  composite  fiber  utilizing  the  following  formula: 

„  _  .  LocalizedHighStress  cr 

btressLoncentration  = - = - 

UniformStressoverRemainder  cravg 


1.  Effect  of  Elastic  Modulus  of  Collagen  on  Stress 
Concentration  of  Composite  Fiber 

Figure  9  shows  how  the  stress  concentration  over  the 
composite  fiber  decreases  as  the  stiffness  of  collagen 
increases.  This  decrease  is  as  expected  because  as  the 
stiffness  of  collagen  increases  its  value  gets  closer  to 
the  stiffness  of  apatite.  From  theory,  the  closer  the 
stiffness  of  the  two  materials,  the  lower  the  stress 
concentration  should  be  (Ugural,  1995)  .  The  results  shown 
in  Figure  9  agree  with  theory. 
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Stress  Concentration  vs  E  collagen 


E  apatite  =  700  GPa  —  —  E  apatite  =  650  GPa  -  -  -  E  apatite  =  600  GPa  E  apatite  =  550  GPa 
—  -  E  apatite  =  500  GPa - E  apatite  =  450  GPa 


Figure  9.  Effect  of  Varying  Elastic  Modulus  of  Collagen 

on  Stress  Concentration 

2.  Effect  of  Elastic  Modulus  of  Hydroxyapatite  on 
Stress  Concentration  of  Composite  Fiber 

Figure  10  shows  that  as  the  stiffness  of  apatite  increases, 
the  stress  concentration  increases  or  maintains  its 
approximate  value.  As  with  Figure  9,  these  results  agree 
with  theory  because  as  you  increase  the  stiffness  of 
apatite  you  are  moving  farther  away  from  the  stiffness  of 
collagen,  thereby  increasing  the  stress  concentration. 
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Stress  Concentration  vs  E  apatite 


E  apatite  (GPa) 


|  -E  collagen  =  10  GPa  —  —  E  collagen  =  15  GPa  -  -  -  E  collagen  =  20  GPa  —  -  E  collagen  =  30  GPa  E  collagen  =  25  GPa~| 

Figure  10.  Effect  of  varying  Elastic  Modulus  of 
Hydroxyapatite  on  Stress  Concentration 


When  analyzing  Figure  9  and  10  together,  we  can  look 
at  the  ratio  of  change  of  the  stress  concentration.  For 
example,  beginning  with  Figure  9  and  a  collagen  stiffness 
of  10  GPa  (with  Eapatite  =  700  GPa)  the  stress  concentration 
is  10.5,  a  50%  change  in  collagen  stiffness  would  give  us 
15  GPa  and  a  stress  concentration  of  9.5.  So  there  is  one 
unit  of  change  in  stress  concentration  with  a  50%  change  in 
collagen  stiffness.  Now,  looking  at  Figure  10,  and  an 
apatite  stiffness  of  450  (Ecollagen  =  15  GPa)  the  stress 

concentration  is  8,  a  50%  change  in  apatite  stiffness  would 
give  us  675  GPa  and  a  stress  concentration  of  9.  Here 
there  is  also  one  unit  of  change  in  stress  concentration 
with  a  50%  change  in  apatite  stiffness. 
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B.  EFFECT  OF  ELASTIC  MODULUS  OF  COLLAGEN  AND 

HYDROXYAPATITE  ON  COMPOSITE  ELASTIC  MODULUS 

Figures  11  and  12  show  the  effect  of  varying  the 

elastic  modulus  of  hydroxyapatite  and  collagen  on  the 
resulting  elastic  modulus  of  the  composite.  As  discussed 
below,  the  changes  in  collagen  have  a  much  greater  impact 
on  the  composite  than  do  the  changes  in  hydroxyapatite. 

1.  Effect  of  Elastic  Modulus  of  Collagen  on  Elastic 

Modulus  of  Composite 

Figure  11  shows  as  the  elastic  modulus  of 

hydroxyapatite  is  kept  constant,  the  elastic  modulus  of  the 
composite  increases  linearly  as  the  elastic  modulus  of 
collagen  increases. 


E  collagen  vs  E  compsite  fiber  with  constant  E  apatite 


- E  apatite  =  700  GPa  —  —  E  apatite  =  650  GPa  -  -  -  E  apatite  =  600  GPa  —  -  E  apatite  =  550  GPa 

E  apatite  =  500  GPa - E  apatite  =  450  GPa 


Figure  11.  Effect  of  Collagen  Elastic  Modulus  on 

Composite  Elastic  Modulus 
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Six  different  values  for  the  hydroxyapatite  elastic  modulus 
are  shown  with  each  line  on  the  graph. 

2.  Effect  of  Elastic  Modulus  of  Hydroyapatite  on 

Elastic  Modulus  of  Composite 

Figure  12  shows  the  slight  increase  of  the  elastic 
modulus  of  the  composite  as  the  elastic  modulus  of 
hydroxyapatite  is  increased. 


E  apatite  vs  E  composite  fiber  with  constant  E  collagen 


E  apatite  (Gpa) 

E  collagen  =  1 0  GPa  —  -  E  collagen  =  20  GPa  ■  -  ■  E  collagen  =  30  GPa  | 


Figure  12.  Effect  of  Hydroxyapatite  Elastic  Modulus  on 

the  Composite  Elastic  Modulus 


The  plot  shows  three  different  values  for  the 
modulus  of  collagen  with  each  different  line.  We 
varying  the  elastic  modulus  of  hydroxyapatite  does 
a  dramatic  effect  on  the  resulting  properties 
composite  material. 


elastic 
see  that 
not  have 
of  the 
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C.  FINITE  ELEMENT  MODEL  COMPARED  TO  AXIAL  DEFORMATION 

CALCULATIONS 

After  utilizing  the  finite  element  model  to  obtain 
values  for  the  elastic  modulus  of  the  composite,  we  wanted 
to  validate  our  simple  analytic  model  using  this  data. 
Figures  13-18  plot  the  FEM  data  with  our  analytic  data  from 
axial  deformation  calculations.  Both  results  compare  very 
well . 

1 .  Constant  Collagen  Elastic  Modulus 

Figures  13-15  show  the  proximity  of  our  finite  element 
model  calculations  to  the  analytic  data  obtained  from  axial 
deformation  calculations  when  varying  the  elastic  modulus 
of  hydroxyapatite.  Errors  are  likely  due  to  approximations 
and  hand  calculations  used  to  obtain  average  stress  results 
from  the  finite  element  model. 

Figure  13,  specifically  shows  the  difference  between 
the  FEM  and  analytic  data  for  Ecoiiagen  =  10  GPa.  Note  that 
the  lines  have  the  same  slope  and  are  within  reasonable 
tolerances  of  each  other. 
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FEM  vs  Analytic  (E  collagen  =  10  GPa) 


E  apatite  (GPa) 


|  —  —  Analytic  FEM  | 

Figure  13.  FEM  vs  Analytic  Data  for  Ecoiiagen  =  10  GPa 

Figure  14,  specifically  shows  the  difference  between 
the  FEM  and  analytic  data  for  Ecollagen  =  20  GPa.  The  lines 
are  almost  right  on  top  of  each  other,  showing  excellent 
correlation  between  the  FEM  data  and  the  analytic 
calculations . 
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FEM  vs  Analytic  (E  collagen  =  20  GPa) 


E  apatite  (GPa) 


—  — Analytic  FEM  | 

igure  14.  FEM  vs  Analytic  Data  for  Econagen  =  20  GPa 

Figure  15,  specifically  shows  the  difference  between 
the  FEM  and  analytic  data  for  Ecouagen  =  30  GPa.  This  graph 
shows  a  bit  more  deviance  from  the  FEM  data  to  the  analytic 
data,  most  likely  due  to  the  error  in  approximating  the 
average  stress  used  to  find  the  composite  elastic  modulus, 
and  due  to  round  off  error. 


FEM  vs  Analytic  (E  collagen  =  30  GPa) 


|  —  —  Analytic  FEM  | 

Figure  15.  FEM  vs  Analytic  Data  for  Ecoiiagen  =  30  GPa 

Overall,  we  determined  from  Figures  13-15,  that  when 
the  elastic  modulus  of  collagen  is  held  constant,  the  FEM 
and  analytic  data  correlate  relatively  well. 

2.  Constant  Hydroxyapatite  Elastic  Modulus 

Figures  16-18  show  the  proximity  of  our  finite  element 
model  calculations  to  the  analytic  data  obtained  from  axial 
deformation  calculations  when  varying  the  elastic  modulus 
of  collagen.  Errors  are  likely  due  to  approximations  and 
hand  calculations  used  to  obtain  average  stress  results 
from  the  finite  element  model. 

Figure  16,  specifically  shows  the  difference  between 
the  FEM  and  analytic  data  for  Eapatite  =  450  GPa.  Again,  we 
see  that  the  values  are  relatively  close. 
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FEM  vs  Analytic  (E  apatite  =  450  GPa) 


Analytic  •  FEM~| 


Figure  16.  FEM  vs  Analytic  Data  for  Eapatite  =  450  GPa 

Figure  17,  specifically  shows  the  difference  between 
the  FEM  and  analytic  data  for  Eapatite  =  550  GPa.  As  with 
Figure  14,  the  lines  are  almost  right  on  top  of  each  other, 
showing  excellent  correlation  between  the  FEM  data  and  the 
analytic  calculat ions . 
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FEM  vs  Analytic  (E  apatite  =  550  GPa) 


Analytic  •  FEM~| 


Figure  17.  FEM  vs  Analytic  Data  for  Eapatite  =  550  GPa 

Figure  18,  specifically  shows  the  difference  between 
the  FEM  and  analytic  data  for  Eapatite  =  650  GPa.  As  with 
Figure  15,  this  graph  shows  a  bit  more  deviance  from  the 
FEM  data  to  the  analytic  data,  most  likely  due  to  the  error 
in  approximating  the  average  stress  used  to  find  the 
composite  elastic  modulus,  and  due  to  round  off  error.  The 
values  are  still  relatively  close. 
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FEM  vs  Analytic  (E  apatite  =  650  GPa) 


Analytic  •  FEM~| 


Figure  18.  FEM  vs  Analytic  Data  for  Eapatite  =  650  GPa 

Overall,  we  determined  from  Figures  16-18,  that  when 

the  elastic  modulus  of  hydroxyapatite  is  held  constant,  the 

FEM  and  analytic  data  correlate  relatively  well. 

D.  EFFECT  OF  VARIOUS  LAMELLA  ANGLES  ON  ELASTIC  MODULUS  OF 
BONE 

Figures  19-24  show  the  elastic  modulus  of  bone  when 
each  layer  of  the  lamellae  is  oriented  at  various  angles. 
For  example,  there  are  12  concentric  lamella,  and  the  first 
lamella  is  oriented  at  15°,  the  second  at  -15°,  the  third  at 
15°,  the  fourth  at  -15°,  and  so  on,  until  you  reach  the 
twelfth  lamella.  Four  possible  orientations  are  shown  in 
the  graphs  in  this  section. 
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1. 


Constant  Collagen  Elastic  Modulus 


Figures  19-21  show  the  elastic  modulus  of  bone  when 
the  lamellae  are  oriented  at  various  angles  and  the  elastic 
modulus  of  collagen  is  kept  constant.  With  all  three 
figures  it  is  obvious  that  the  overwhelming  contributor  to 
the  elastic  modulus  of  bone,  is  the  elastic  modulus  of 
collagen . 


E  bone  effected  by  lamella  angles  (E  collagen  =  10  GPa) 


—  700/10  - 650/10  -  -  -600/10  550/10  —  -  500/10  450/10 


Figure  19.  Bone  Elastic  Modulus  with  Ecoiiagen  =  10  GPa 

Figure  19,  specifically  addresses  the  case  when  Ecoiiagen 
=  10  GPa.  We  see  that  the  various  values  for  the  elastic 
modulus  of  bone  are  similar  for  each  value  of 
hydroxyapatite,  but  as  the  degree  of  orientation  is 
increased,  the  bone  elastic  modulus  decreases,  as  expected. 
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The  increasing  angle  of  orientation  moves  each  lamella 
further  away  from  its  strictly  longitudinal  configuration, 
and  therefore  causes  the  bone  elastic  modulus  to  continue 
to  decrease.  With  collagen  being  so  soft,  the  stiffness  of 
bone  remains  relatively  low  for  each  value  of  Eapatite- 


E  bone  effected  by  lamella  angles  (E  collagen  =  20  GPa) 


700/20  - 650/20  -  -  -600/20  550/20  —  -  500/20  450/20 


Figure  20.  Bone  Elastic  Modulus  with  Ecoiiagen  =  20  GPa 

Figure  20,  specifically  addresses  the  case  when  Ecoiiagen 
=  20  GPa.  The  same  trend  in  Ebone  is  recognized.  The 
further  away  from  zero  degrees  the  lamella  orientation 
gets,  the  lower  the  value  for  bone  elastic  modulus.  We  can 
see  what  a  difference  the  harder  collagen  (E  =  20  GPa  vice 
E  =  10  GPa)  makes  on  the  bone  elastic  modulus  when 
comparing  Figure  20  to  Figure  19. 
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E  bone  effected  by  lamella  angles  (E  collagen  =  30  GPa) 


Figure  21.  Bone  Elastic  Modulus  with  Ecoiiagen  =  30  GPa 

Figure  21,  shows  results  with  Econagen  =  30  GPa.  When 
raising  collagen's  elastic  modulus  this  high,  we  see  a 
marked  increase  in  the  slope  of  each  line,  and  therefore  a 
greater  increase  in  the  property  of  bone.  Although,  when 

getting  further  away  from  the  zero  degree  orientation 

(i.e.,  the  60°  orientation)  the  values  for  the  stiffness  of 
bone  are  very  similar  to  previous  values  with  the  softer 
collagen.  This  is  due  to  the  high  degree  of  orientation 

for  each  concentric  lamella. 

2 .  Constant  Hydroxyapatite  Elastic  Modulus 

Figures  22-24  show  the  elastic  modulus  of  bone  when 
the  lamellae  are  oriented  at  various  angles  and  the  elastic 
modulus  of  hydroxyapatite  is  kept  constant.  With  all  three 
figures  we  see  the  same  trend  of  bone  stiffness,  and 
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realize  that  the  major  contributor  to  its  properties  are 
that  of  collagen.  Note  the  striking  similarities  between 
figures  22,  23,  and  24. 


E  bone  effected  by  lamella  angles  (E  apatite  =  450  Gpa) 


Figure  22.  Bone  Elastic  Modulus  with  Eapatite  =  450  GPa 

Figure  22  shows  the  specific  case  of  Eapatite  =  450  GPa. 
As  in  the  previous  figures  we  see  a  decreasing  trend  of  the 
elastic  modulus  of  bone  as  the  degree  of  orientation  gets 
greater  and  greater.  This  is  again  due  to  the  greater 
distance  from  the  longitudinal  orientation. 
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E  bone  effected  by  lamella  angles  (E  apatite  =  550  Gpa) 


Figure  23.  Bone  Elastic  Modulus  with  Eapatite  =  550  GPa 

Figure  23  shows  the  specific  case  of  Eapatite  =  550  GPa. 
Remarkably  similar  to  Figure  22,  providing  an  example  of 
how  the  properties  of  bone  are  not  reliant  on  the 
properties  of  hydroxyapatite,  but  instead  on  the  properties 
of  collagen. 

Figure  24  depicts  the  same  trend  in  results,  but  with 
Eapatite  =  650  GPa . 
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E  bone  effected  by  lamella  angles  (E  apatite  =  650  GPa) 


—  -  650/10  —  -  650/15  -  -  -  650/20  650/25  - 650/30 

Figure  24.  Bone  Elastic  Modulus  with  Eapatite  =  650  GPa 
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V. 


CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  objective  of  this  study  was  to  derive  an  analytic 
model  of  the  basic  hierarchical  structure  of  the  human 
bone.  The  model  computed  the  stiffness  of  composite 
collagen  fibers  comprised  of  collagen  fibrils  and 
hydroxyapatite  mineral  crystals.  Next,  the  stiffness  of  the 
concentric  lamella  was  computed  utilizing  the  stiffness  of 
the  collagen  fibers  and  layer  information.  Finally,  the 
effective  stiffness  of  the  bone  was  estimated.  In  order  to 
determine  the  stiffness  of  the  collagen  fiber  and  validate 
the  analytic  model,  a  three-dimensional  finite  element 
model  was  developed.  Once  the  model  was  validated,  a 
series  of  parametric  studies  were  conducted  to  understand 
what  parameter  (s)  affected  the  stiffness  of  the  bone  most 
significantly . 

The  results  of  the  parametric  study  demonstrate  that 
the  elastic  modulus  of  collagen  is  much  more  vital  in 
determining  the  stiffness  of  bone  than  is  the  elastic 
modulus  of  hydroxyapatite.  Both  materials  had  a  relatively 
similar  effect  on  the  stress  concentration  of  the  composite 
fiber.  Therefore,  if  one  were  to  construct  an  artificial 
human  bone  they  would  need  to  pay  close  attention  to  the 
material  used  to  replicate  collagen,  and  should  not  be  so 
concerned  about  the  material  used  for  hydroxyapatite. 
Another  important  result,  was  that  the  lower  the  degree  of 
orientation  from  the  centerline,  the  higher  the  elastic 
modulus  of  bone.  But  what  was  not  analyzed  in  this  study 
was  the  effect  of  torsion  on  the  bone  stiffness.  To  combat 
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the  effects  of  torsion,  would  require  a  higher  degree  of 
orientation  of  the  lamella. 

B .  RECOMMENDATIONS 

The  following  recommendations  for  future  research 
include : 

The  current  finite  element  model  was  analyzed  for  the 
axial  load.  Future  work  could  consist  of  analyzing  the 
model  with  different  loading  conditions  including  a 
torsional  and  a  bending  load. 

This  research  utilized  the  elastic  modulus  of 
hydroxyapatite  and  collagen  in  determining  the  elastic 
modulus  of  bone.  In  the  future,  other  bone  properties  such 
as  the  transverse  modulus  and  shear  modulus  could  be 
analyzed  to  determine  if  the  same  trends  occur.  The 
interface  strength  between  the  collagen  and  apatite  could 
be  another  important  parameter  to  analyze. 
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APPENDIX  A  -  LONGITUDINAL  ELASTIC  MODULUS  OF 
COMPOSITE  FIBER  MATLAB  PROGRAM 


%This  program  is  used  to  calculate  the  total  stiffness 
%coefficient  for  one  section  of  collagen  fiber  (with 
%apatite)  as  well  as  the  Elastic  Modulus  of  the  composite 
%f iber . 


%  kl  k2  k3 

%Theory 

%1/k  =  1 /kl  +  l/k2  +  l/k3 

%1/k  =  (k2k3  +  klk3  +  klk2 ) /  (klk2k3 ) 


clear 

clc 

r  =  60e-9;  %radius  in  meters 

Ec  =  10e9;  %Modulus  of  collagen  in  Pascals  (N/m^2) 

LI  =  8e-9;  %length  of  initial  and  final  section  of 

collagen  in  meters 

L2  =  48e-9;  %length  of  apatite  in  meters 

n  =  12;  %nr  of  pieces  of  apatite  around 

circumference 

Aa  =  4 . 97955e-9*4*2e-9;  %Area  of  apatite  crystal  (edge 

length* 4* width) 

Ea  =  700e9;  %Modulus  of  apatite  in  Pascals  (N/m^2) 

L  =  2*L1  +  L2;  %Total  length  of  composite 

%kl  =  (AE/L1 ) 

kl  =  (pi*r/'2*Ec)  /LI; 

alpha  =  0.25; 

%k2  =  k21  +  k22  (collagen  +  apatite) 
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%k2  =  (AE/L)c  +  (AE/L) a 
k2a  =  alpha* (n*Aa*Ea) /L2 ; 
k2c  =  (pi*r/'2*Ec)  /L2; 
k2  =  k2a  +  k2c; 
k3  =  kl; 

k  =  (kl*k2*k3) /  (  (kl*k2)  +  (k2*k3)  +  (kl*k3)); 

%k  =  (AE/L) total 
Ecomposite  =  (  (k*L)  /  (pi*r/'2 )  ) 
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APPENDIX  B  -  ELASTIC  MODULUS  OF  COMPOSITE  FIBER 
ORIENTED  AT  DIFFERENT  ANGLES  MATLAB  PROGRAM 


%To  determine  modulus  of  elasticity  for  collagen  at  15,  30, 

%45,  60,  -15,  -30,  -45,  -60  degree  angles. 

clc 

clear 

nu  =  0.3; 

Efiber  =  33e9; 
data  =  [ ] ; 

1=0; 

%Try  angles  of  15,30,45,60,-15,-30,-45,-60  degrees 

for  theta  =  [(pi/12)  (pi/6)  (pi/4)  (pi/3)  (-pi/12)  (-pi/6) 

(-pi/4)  (-pi/3) ] ; 


Qll  =  Efiber/  (l-nuA2)  ; 

Qllbar  =  Qll*  (cos (theta) ) A4; 

Q12bar  =  Qll*  (sin  (theta)  )  A2*  (cos  (theta)  )  A2; 
Q22bar  =  Qll*  (sin (theta) ) A4; 

Q16bar  =  Qll*  (sin  (theta) )* (cos (theta) ) A3; 
Q26bar  =  Qll* (sin (theta) ) A3* (cos (theta) ) ; 
Q66bar  =  Qll* (sin (theta) ) A2* (cos (theta) ) A2; 


%Qbar  =  [Qllbar  Q12bar  Q16bar;  Q12bar  Q22bar  Q26bar;  Q16bar 
Q2 6bar  Q66bar ] ; 

Qbar  =  zeros  (3, 3);  %Set  matrix  of  zeros  as  most 
%values  are  assumed  zero 


Qbar  (1,1)  =  Qllbar; 
Qbar(2,2)  =  Q22bar; 


Force  value  for  (1,1) 
Force  value  for  (2,2) 


theta 


Qbar; 
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Sll  = 

1/Qbar (1, 1) ; 

S22  = 

1/Qbar  (2,2) ; 

Elbar 

=  1/S11 

%Value  for 

El 

E2bar 

=  1/S22; 

%Value  for 

E2 

i  =  i+1; 

data  =  [data;  Qllbar  Q22bar] ; 
disp  (  [ '  Qllbar  Q22bar']); 
data; 


end 
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APPENDIX  C 


-  LAMELLA  AND  BONE  ELASTIC  MODULUS 
MATLAB  PROGRAM 


%Assuming  that  there  are  12  layers  of  collagen  which  make 
%up  the  lamella  and  that  each  layer  is  120  nm  thick.  Take 
%the  values  for  E  from  the  "transf ormedmatrix"  m-file  for 
%various  angles. 


clc 

clear 


%theta  =  pi/12  (15  degrees) 

Elbarl5  =  31.568e+09; 

Elbarnl5  =  Elbarl5; 

%theta  =  pi/6  (30  degrees) 

Elbar30  =  20.398e+09; 

Elbarn30  =  Elbar30; 

%theta  =  pi/4  (45  degrees) 

Elbar45  =  9.0659e+009; 

Elbarn45  =  Elbar45; 

%theta  =  pi/3  (60  degrees) 

Elbar60  =  2.2665e+009; 

Elbarn60  =  Elbar60; 

%radius  calculations  for  each  of  12  layers  (assuming  radius 
%of  each  layer  is  equivalent  to  diameter  of  a  single 
%collagen  fiber) 

rl  =  60e-9; 

r2  =  rl  +  120e-9; 

r3  =  r2  +  120e-9; 
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r4  = 

r3 

+ 

120e-9; 

r5  = 

r4 

+ 

120e-9; 

r6  = 

r5 

4 

12  Oe- 9; 

r7  = 

r  6 

4 

120e-9; 

r8  = 

r7 

4 

12  Oe- 9; 

r  9  = 

r  8 

+ 

120e-9; 

rl  0 

=  r  9 

1  4 

120e-9; 

rl  1 

=  r  1 0 

4  120e-9; 

rl2 

=  r  1 1 

4  120e-9; 

%Area  calculations  for  each  layer 

A1  =  pi*rlA2; 

A2  =  pi*  ( r2  A2-r 1 A2  )  ; 

A3  =  pi*  ( r3  A2-r2  A2 ) ; 

A4  =  pi*  ( r4  A2-r3  A2  )  ; 

A5  =  pi*  ( r5  A2-r4  A2  )  ; 

A6  =  pi* (r6A2-r5A2) ; 

A 7  =  pi*  ( r7  A2-r 6 A2 ) ; 

A8  =  pi*  ( r 8 A2-r7 A2 ) ; 

A9  =  pi* ( r 9 A2-r 8  A2 ) ; 

A10  =  pi* ( r 1 0  A2-r 9 A2 ) ; 

All  =  pi*  (rllA2-rlOA2) ; 

A12  =  pi* (rl2A2-rllA2) ; 


%Total  area  of  lamella 
Atotal  =  pi*rl2A2; 

Atotalsum  =  A1+A2+A3+A4+A5+A6+A7+A8+A9+A1 0+A1 1+A12 

%Total  Modulus  of  Elasticity  (El) 

%Try  opposite  layers  (15,-15,15,-15...) 
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Eltotall5= (Al*Elbarl5) + (A2*Elbarnl5) + (A3*Elbarl5) + (A4*Elbar 
nl5) + (A5*Elbarl5) + (A6*Elbarnl5) + (A7*Elbarl5) + (A8*Elbarnl5) + 
(A9*Elbarl5) + (A10*Elbarnl5) + (All*Elbarl5) + (A12*Elbarnl5) ) /A 
total ; 


%Multiply  by  0.7  to  account  for  70%  of  bone 
Elbonel5  =  Eltotall5*.7 


%Try  opposite  layers  (30,-30,30,-30...) 

Eltotal30= ( (Al*Elbar30) + (A2*Elbarn30 ) + (A3*Elbar30) + (A4*Elba 
rn30) + (A5*Elbar30) + (A6*Elbarn30 ) + (A7*Elbar30) + (A8 *Elbarn30 ) 
+ (A9*Elbar30) + (A10*Elbarn30) + (A1 1 *Elbar30 ) + (A12 *Elbarn30 ) ) / 
Atotal ; 


%Multiply  by  0.7  to  account  for  70%  of  bone 
Elbone30  =  Eltotal30*.7 


%Try  opposite  layers  (45,-45,45,-45...) 

Eltotal45= (Al*Elbar45) + (A2*Elbarn45) + (A3*Elbar45) + (A4*Elbar 
n45) + (A5*Elbar45) + (A6*Elbarn45) + (A7*Elbar45) + (A8*Elbarn45) + 
(A9*Elbar45) + (A10*Elbarn45) + (All*Elbar45) + (A12*Elbarn45) ) /A 
total ; 


%Multiply  by  0.7  to  account  for  70%  of  bone 
Elbone45  =  Eltotal45*.7 


%Try  opposite  layers  (60,-60,60,-60...) 

El total  60= ( (A1 *Elbar 60 )  +  (A2 *Elbarn60 )  +  (A3*Elbar60)  +  (A4*Elba 
rn60) + (A5*Elbar60) + (A6*Elbarn60 ) + (A7*Elbar60) + (A8 *Elbarn60 ) 
+ (A9*Elbar 60 ) + (A1 0 *Elbarn60 ) + (A1 1 *Elbar 60 ) + (A12 *Elbarn60 ) ) / 
Atotal ; 


%Multiply  by  0.7  to  account  for  70%  of  bone 
Elbone60  =  Eltotal60*.7 
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%Try  random  layers  (15,-15,30,-30,45,-45,15,-15,30,-30,45,- 
%4  5 ) 

Eltotalrl= ( (Al*Elbarl5) + (A2*Elbarnl5) + (A3*Elbar30) + (A4*Elba 
rn30) + (A5*Elbar45) + (A6*Elbarn45) + (A7*Elbarl5) + (A8*Elbarnl5) 

+ (A9*Elbar30) + (A10*Elbarn30) + (All*Elbar45) + (A12*Elbarn45) ) / 
Atotal ; 


%Multiply  by  0.7  to  account  for  70%  of  bone 
Elbonerl  =  Eltotalrl*.7 


%Try  random  layers  (45,-45,  60,-60,45,-45,  60,-60,45,-45,  60,- 
%60 ) 

Eltotalr2= ( (Al*Elbar45) + (A2*Elbarn45) + (A3*Elbar60) + (A4*Elba 
rn60) + (A5*Elbar45) + (A6*Elbarn45) + (A7*Elbar60) + (A8 *Elbarn60 ) 

+ (A9*Elbar45) + (A10*Elbarn45) + (A1 1 *Elbar 60 ) + (A12 *Elbarn60 ) ) / 
Atotal ; 


%Multiply  by  0.7  to  account  for  70%  of  bone 
Elboner2  =  Eltotalr2*.7 
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